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that the behavior is quite different if the node always enters in the dynamics as its own input 
(self-regulation) or not. The same conclusion holds for the Kauffman NK model. Moreover, the 
simulation results and the mean-field ones (i) agree well when there is no self-regulation, and (ii) 
disagree for small p when self-regulation is present in the model. 

PACS numbers: 05.10.-a, 05.45.-a, 87.18.Sn 
Keywords: 



I. INTRODUCTION 

Some physicists have claimed that it is possible to 
roughly classify its science branches in physics of small, 
big and finally complex systems. Although such clas- 
sification may appear simple, it suits very well for the 
present work. Complex systems are known as entities 
composed of a large number of agents sharing a rich 
set of simple and non linear interactions. In such sys- 
tems, different behaviors can be achieved if the interac- 
tion change, even if the interacting agent remains the 
same. These systems are well modeled by using network 
concepts, where the agents are called nodes, and interac- 
tions appear as links. 

Since the release of the Barabasi and Albert paper |l| 
about growing networks, a lot of knowledge has been 
achieved regarding the topology and properties of such 
objects In addition, it was found that these types 

of networks can be found in a variety of fields like social 
relation, voting, disease spreading, the WWW, neural 
and regulatory networks [3-ll5|. 

The interaction among nodes of a network can be mod- 
eled in very different ways. In the most simple models, 
it is assumed that (i) nodes can only have two states 
(off) or 1 (on), and that (b) the dynamics is done via 
Boolean functions. These models, very suitable for sim- 
ulations in silico, can be be very useful to study agents 
that interact using on/off states as happens in gene ex- 
pression/repression 10\ and protein activation/inhibition 
|2j. Models involving Boolean dynamics are called Ran- 
dom Boolean Networks (RBNs). The first RBN model 
was introduced by Stuart Kauffman, and it is known as 
NK model since it was composed by N nodes, each one 
with its own random K inputs, and with the Boolean 



functions randomly chosen |T7I |. 

In general, all RBNs share some common features. 
They are directed networks with N nodes, and a node i 
has ki inputs that regulate its next state via some random 
Boolean function Gi. The topology of the interactions 
can be expressed by using a proper connectivity distri- 
bution of inputs P(k) and the type of the Boolean func- 
tions can be changed in order to adequate the restrictions 
imposed by the problem. The Boolean functions can be 
chosen basically in two different ways: (a) all functions 
are chosen a priori, and they are kept constant during 
time evolution - the so-called quenched models; and (b) 
the functions change during the dynamics, meaning that 
for each time step, each node have a new function - the 
so-called annealed models. Another important dynami- 
cal feature is the update of the nodes. In the synchronous 
or parallel mode, all nodes are updated simultaneously in 
each time step. On the other hand, in the asynchronous 
or serial mode, first we randomly choose a node that 
is updated immediately; then this procedure is repeated 
until we have updated N nodes in a single time step. 

Recently scale-free distributions came up in the arena 
of RBNs in order to match biological network scenar- 
ios. The first works have only used computer simulations 
P^ - |20| . Later, some authors have focused on analytical 
approaches |2lT[24| giving a more detailed insight of the 
problem. In this work we study random Boolean dy- 
namics in scale-free networks. In our model we have 
nodes with Boolean variables (0,1). The connections 
among nodes obey the topological structure of scale-free 
networks and the relation among nodes variables is per- 
formed via randomly chosen Boolean functions. We sup- 
pose that the dynamics is driven only by XOR and AND 
functions, which appear for each node with probabili- 



2 



ties p and 1 — p, respectively. We chose the AND and 
XOR logical functions in order to simplify the model, 
avoiding the necessity of defining 2 2 different Boolean 
relations for each node. Note that any Boolean function 
can be written as a linear combination of AND, XOR 
and OR. The chosen functions are examples of extreme 
cases. When the AND function is applied to a set with 
K Boolean variables, we will obtain 1 only if all variables 
are l's. On the other hand, when the XOR function is 
applied to the same set, we will obtain 1 if the number 
of l's of the set is odd. These observations imply that 
AND represents a very selective dynamical rule, only one 
configuration of the 2 K possible ones furnishes 1 as out- 
put, while XOR is related to a non-selective dynamical 
one because half of the set configurations give 1 as out- 
put. Moreover, it is known from previous works that 
the AND function leads to an ordered regime with two 
fixed points, where all variables are O's or l's, and that 
the dynamical behavior generated by the XOR function 
is more complex. Since this kind of functions are the 
Boolean counterparts of real reactions in cell regulatory 
system [2^, [2(| , they are important to the study of bio- 
logic networks. We choose the parallel mode as updated 
method. In order to compare with another models, and 
for technical reasons, we also study such dynamics in net- 
works without a scale-free topology. The Kauffman NK 
model is briefly discussed as well. We performed two dis- 
tinct types of dynamics. In the first case the node that 
will be updated is regulated only by the nodes which 
are connected to it. Rarely the node is connected to it- 
self. It means that its state is defined by the state of its 
neighbors, and its own state is almost never taken into 
account. In the second case we explore self-regulation, 
which means that the new state of the updating node 
is defined by its neighbors and always for its own state. 
We choose to study such case because self-regulation is a 
well know feature of genetic regulatory networks [27| . In 
section II we introduce the scale-free networks and the 
Boolean dynamics used in this work. The numerical sim- 
ulations of the scale- free networks are discussed in section 
III. In section IV we present a mean-field (annealed) ap- 
proximation for these dynamics, leading to an analytical 
way to calculate the average density of l's and the Ham- 
ming distance. We apply the annealed approximation to 
networks without a scale-free topology, to the Kauffman 
NK model and to scale-free networks. The comparison 
between the numerical results and those of the annealed 
approximation is presented in section V. We summarize 
our results in the last section. 



II. NETWORKS AND THE DYNAMICS 

We have generated two classes of distinct networks, 
classified by the smallest number of links that a node 
can have, in other words, by the smallest possible connec- 
tivity k m i n of the network. These networks were grown 
by the Growing Network with Re-direction algorithm [28j 



and can be classified as 

• k m i n = 1 - a node of the network, called old node, 
is selected with uniform probability; then a new 
node is linked to it with probability 1 — r or it 
is redirected to the ancestor of the old node with 
probability r; 

• k m i n = 2 - a new node has two links; the first link 
with an old node, selected with uniform probabil- 
ity, is established with probability 1 — r or it is redi- 
rected to one of the two ancestors of the old node 
with probability r; the same procedure is repeated 
for the second link. 

For k m i n = 1, initially we have three nodes cyclically 
connected (the ancestors of nodes 1, 2 and 3 are 3, 1 and 
2 respectively) . A new node is randomly connected to an 
old node (one of the three initial nodes) . Then, this new 
link can be redirected to the ancestor of the old node with 
probability r. This growing algorithm is repeated until 
we have N nodes in the network. When k m i n = 2, each 
one of the initial three nodes has the other two nodes 
as ancestors. Now, a new node is randomly connected 
to two old nodes, and each new link can be redirected to 
the ancestor with probability r. We repeat this procedure 
until we have a network with N nodes. These algorithms 
create scale- free networks characterized by a connectivity 
distribution P{k) ~ fc~ 7 with 7 = r" 1 + 1 When 
r = 0, the nodes are linked in a entirely random fashion 
and for r = 1, all nodes are connected to one of the 
three initial nodes (super- hubs). The linear preferential 
attachment model of Barabasi and Albert is obtained 
for r = 0.5 (7 = 3). In this work, we consider three 
typical values of r: r = 0.5 (the Barabasi and Albert 
model), r — 0.8, representing models with large hubs 
and r = 0.35, representing models with small hubs. 

A logical variable o~i(t) is assigned to each node i and 
the state of the network at time t is represented by a set 
of Boolean variables (<7i (t) , a% (t) , 03 (t) , crjv(i))- Each 
variable o~i{t) is controlled by fej elements of the network 
W(t)}ki = {cr il (i), cr i2 (t), cr lfc . (t)}. If k min = 2, ki is 
the connectivity of z-th node and the control elements 
are the nodes connected to it. When k m - m = 1, we have 
nodes with only one link. Since we need two inputs to 
apply the Boolean functions, it is natural to assume that 
the node itself must always participate of the dynamics. 
Now, the control elements are the z-th node itself and 
nodes connected to it. It means that each node has an 
extra link to itself (self-regulation) . In this case, ki is the 
connectivity plus 1. The dynamics is given by 

a t (t + l) = G l ({a(t)} ki ), (1) 

where Gi is the random function 

G _ I AND with probability 1-p „ 
[ XOR with probability p, 
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that is assigned to each node i. Here p is an external 
parameter that controls how the logical functions AND 
and XOR are distributed in the network. 

An initial state {c(0)} is created by assigning ran- 
domly O's and l's to all nodes. A damaged copy {cr(0)} 
of the initial state is also created by changing the value of 
only one randomly chosen node. Since the Hamming dis- 
tance of two configurations is the number of nodes that 
have different values in these conf iguratio ns, the Ham- 
ming distance between {cr(0)} and {c(0)} is 1. Both the 
initial state and its copy evolve under the control of equa- 
tion ([T]) . Once the new state of all nodes is calculated the 
entire network is updated (synchronous update) and the 
system goes to the next Monte Carlo time step (mcs). 

We characterize the dynamical behavior by the average 
density of l's 

M(p,t)=jta |>(*>). 

and by the average of the Hamming distance 

N 

Y v \ .V 



i=l 



Here, (. . .) is an average over the initial states ({cr(0)}, 
{c(0)}), and over sets of links of a grown network with 
a specific 7 and with the same p. 

After a transient time these quantities reach the sta- 
tionary values M(j>) and D{p) that can be defined as 



M(p) 



D{p) 



lim 



lim 



t+T 



M(p,t')dt', 



t+T 



D(p,t')dt'. 



(3) 



(4) 
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FIG. 1: Log-log plots of M(p) vs p for: (a) a network with 
kmin = 1, r = 0.5 and different N (b) a network with k m i n = 
1, N = 10 4 and different r 




III. NUMERICAL SIMULATIONS 



FIG. 2: Log-log plot of M(p) vs p for k m in = 1, N varying 
from 1 x 10 4 up to 4 x 10 4 and r = 0.35, 0.5 and 0.8 . It shows 
the best fit of the coalesced sets. 



A. Data and results for k min — 1 

In order to consider finite size effects, we grew networks 
with N = 1 x 10 4 , N = 2 x 10 4 and N = 4 x 10 4 nodes. 
The averages were performed with a number of samples 
varying from 10 2 (large p and large N) up to 5 x 10 4 
(small p and small N). The probability p was taken in 
the interval [0.001, 0.8] for the three values of r (0.35, 
0.5 and 0.8). Since for p < 0.001 the averages quantities 
were very small, we decided that a good lower limit was 
p = 0.001. The upper limit p — 0.8 was chosen because 
the values of the average quantities were similar, in a log- 
log scale. The stationary values, M(jp) and D(p), were 
reached after a very short transient time (20 mcs). Esti- 
mations of M(p) and D(p), which are defined in Eqs. Q 
and were performed by considering t = 20 mcs and 



T = 80 mcs. These stationary values remain basically 
the same if we increase both t and T. 

We can see from Fig. Q] that M(p) behaves as a power 
function of type p m for the entire range of p. By compar- 
ing the behavior of the smallest network (N = 1 x 10 4 ) 
with the largest one ( N — 4 x 10 4 ) we can see that fi- 
nite size effects are small for M{p). Moreover, it seems 
that the exponent m does not depend on r. In order to 
evaluate the exponent m we coalesce all different sets [N 
varying from 10 4 up to 4 x 10 4 and r from 0.35 up to 0.8) 
and we do a best fit. This is shown in Fig. [5J We obtain 
that 

M(p) = ap m , 

with a = 0.46 ± 0.01 and m = 0.96 ± 0.01. Note that 
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FIG. 3: (a) Log-log plot of D(p) vs p for a network with 
kmin = 1, r = 0.5 and different TV; (b) Linear-log plot of 
D(p) vs p for a network with k m i n = 1, N = 2 x 10 4 and 
different r. 




FIG. 4: (a) Linear-log plot of D(p) + B vs p for a network 
with k m in = 1, r = 0.5 and different N; (b) Log-log plot of 
D(p) + B vs p for a network with k m i n = 1, AT = 2 X 10 4 and 
different r. 



we have evaluated the exponent by considering approxi- 
mately 2 orders of magnitude in the p variable and that 
the fit is very good. In fact, in all fitted data, we obtained 
a correlation coefficient larger that 0.999. 

Plots of D{p) versus p are shown in Fig. [3] for r — 0.5 
and networks with different sizes. We can see that D(p) 
has a power law behavior (D(p) ~ p d ) only when the 
probability p is close to p = 0. Outside of these region, 
D(p) grows exponentially. By comparing the behavior of 
networks with different sizes, we observe that finite size 
effects are now important. 
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FIG. 5: Log-log plot of M(p) vs p for k min = 2, N varying 
from 1 x I0 4 up to 4 x I0 4 and r = 0.2, 0.5 and 0.8 . It shows 
the best fit of the coalesced sets. 



It turns out that all results can be well fitted by 

D{p) = J Bexp(Cp) - B : 

with C and B depending on the size of the network and 
on the parameter r. This can be seen in Fig. 2J When 
p w 0, we have that D ~ BCp. 

B. Data and results for k m m = 2 

The simulation for the networks with k m i n = 2 were 
realized in the same way that for k m i„ = 1, and the 
probability p was taken in the interval [0.01, 0.9] for r = 
0.2, 0.5 and 0.8. The range of p is rather narrow in this 
case, since the dynamics for k m i n = 2 is more sensible to 
p values, leading to M(p,t) — > and D(p,t) — > when 
p ~ 0.001. This behavior perhaps is related to the one 
found for RBN with only part of the canalyzing functions 
as update functions, K = 2, and small p + , where pu- 
is the probability that a connection be excitatory [29j. 
Fig. [S] is similar to fig. [21 where all sets of r and N 
are coalesced in one plot. We can see that the finite size 
effect is very small and we have 

M(p) = ap m , 

with a = 0.70 ± 0.02 and m = 1.79 ± 0.01. 

As we can see in Fig. [SJ the Hamming distance D (p) 
follows a linear dependence with p for large values of p 
(the plot shows p > 0.3). In this region the behavior of 
D(p) is almost independent of r, and we have: 

Dip) ~ ap, 

where a = 0.78 ± 0.02, 0.76 ± 0.02 and 0.72 ± 0.02 for 
r = 0.2, 0.5 and 0.8 respectively. However, the behavior 
of Dip) for p < 0.3 does not provide any suitable fit 
since for low values of p we have a high concentration of 
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FIG. 6: Linear plot of D(p) vs p for kmin 
and 0.8 for N — 10 3 . 



= 2, r = 0.2, 0.5 



AND function and a low concentration of XOR function, 
leading most of nodes to the state. 

Note that results discussed in this section are valid for 
other values of r. We discuss only the cases r = 0.5 
(linear preferential attachment), one case with r < 0.5 
(r = 0.35 or r = 0.2) and one case with r > 05 (r — 
0.8 or r = 0.7) because they represent typical behaviors. 
We have simulated other cases, with less samples, and 
the results are similar. Although we present D(p,t) for 
the initial condition D(p, 0) = 1, we have also simulated 
cases with D(p, 0) > 1. We find similar results, probably 
because this new initial condition is a later state of the 
initial condition with the smallest Hamming distance. In 
the next section we will develop a mean- field approach in 
order to see if its results agree with the numerical ones 
just obtained. 



IV. MEAN-FIELD APPROXIMATION 

In this section we present a mean- field approach (MF). 
It is based on a work of Derrida and Pomeau [301 ] . in 
which an annealed approximation was done for the Kauff- 
man NK model. The NK model is a cellular automaton 
with N nodes holding logical variables. Each node Oi is 
connected with another K nodes of the network meaning 
that all nodes have the same connectivity K. The dy- 
namics is given by Eq. ([T]) where Gj is a random Boolean 
function. Although the model is defined by quenched dis- 
order, i.e, the Boolean function Gi and the K nodes con- 
nected to each node Gi are only randomly chosen at the 
initial time, in the Derrida and Pomeau approximation 
it is assumed an annealed disorder. It means that the 
Boolean function and the K nodes are randomly chosen 
at each time step. Moreover, in such approximation the 
effect of the Boolean functions on a node is described by 
probabilities that the output be and 1. The approxi- 
mation to our problems is similar to that of Derrida and 
Pomeau. However, a difference appears in the application 



of the Boolean functions. Instead of using a probabilis- 
tic description for the effect of the Boolean functions, we 
determine the effect of applying the XOR and AND op- 
erators in each configuration. We will first discuss the 
case k m i n — 2 because it is more illustrative than the 
more simple case fc m in = !■ 



A. Average density of l's for k m i n — 2 

Let us start our evaluation for the model without a 
scale-free topology. The dynamics is given by equations 
(HJ) and @ . Suppose that the configuration of the system 
at time t, {(7i(i)}, consists of n nodes with a = I and 
N — n nodes with a = 0. In order to study the density 
of l's, we can separate the configuration {<Ji(t)} in two 
sets: (i) set A(t) where all the nodes have a — 0, and (ii) 
set Bit) where all the nodes have a = 1. Now we must 
randomly choose K nodes that are linked to the node i. 
The probability that a given link comes from set A{t) is 
1 — x = (N — n)/N and the probability that it comes from 
B(t) is x — n/N. Since each node has K links, the list 
Pk of the probabilities of the possible link configurations 
is 



Pk = 



(l-x) K , 



K 

K -2 



K 
K - 1 

K-2 2 



(l-x) K - z x z , 



(1-x) 

K 




K-l„ 



(5) 



„K 



We must now evaluate the output of each possible con- 
figuration under the application of operators AND and 
XOR. The Boolean function AND generates an output 1 
if all inputs come from set Bit). The probability of this 
configuration is given by the last term of the list (JSJ . The 
function XOR, however, produces 1 as output when its 
number of inputs equal to 1 is an odd number, meaning 
that the number of links coming from set Bit) is odd. 
Therefore the probability X of obtaining 1 as output is 



K 



m—1 

m odd 



X=p^2 ( K )(l-x) K - m x m +qx K , 



where p and q = 1 — p are the probabilities related to the 
XOR and AND operators. Since we have that 



A" 



the sums of the even and odd terms are given by 

K 



J2 [ K ^ y K-m x m iy + X) K + iy - X) K 

m=0 

m even 



rn 



(6) 



E 

m—l 

m odd 



m J 2 
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It follows that X = |[1 - (1 - 2x) K ] + qx K . Assuming 
homogeneity, we can identify X as M(p,t + 1), the frac- 
tion of l's at time t + 1. Then we obtain that M(p, t + 1) 
depends only on M(j>, t) as 

M(p,t+l) = |{1 - [1 - 2M(p,t)] K } + qM(p,tf. (8) 

Finally let us consider the scale-free topology. Each 
node now has k links with probability P(k). Therefore 
the above equation can be written as 



M(p,t + 1) = | jl-^P(fe)[l-2M(p,t)] fc j 

OO 

+ g£p(*)M(p,t) fc . (9) 



fe=2 



B. Average Hamming distance for k m in = 2 

The calculation of the Hamming distance is done in 
a similar way. Let us again first study the model with- 
out scale-free topology. At time t we are interested in 
the configuration {cr(t)} resulting from the evolution of 
the initial configuration, and in the configuration {a(t)}, 
which appears from the evolution of the initial damaged 
copy. Suppose that they differ by n nodes. Following 
Derrida and Pomeau [3fJ], we define two sets: £{t) and 
J-(t). The first one is the set of all nodes of {cr(t)} and 
{cr(t)} that have the same values. The set J-(t) is com- 
posed by the n nodes that have different values in the 
two configurations. Therefore the nodes which have all 
links coming from set £{t) will have the same values at 
time t+1 in the {c(i)} and {a(t)} configurations and 
they will not contribute to the Hamming distance. On 
the other hand, the Hamming distance could be changed 
by the nodes which have at least one link coming from 
the n nodes of .F(i). 

Let us also define Eq and E\ as the number of nodes 
in the set £ (t) with the values and 1, respectively. Fq is 
the number of nodes in the set F{t) that have a(t) = 0, 
and Fi is the number of nodes with a(t) = 1. Observe 
that F + Fi = n and Eq + E\ = N — n. The next step 
is to focus in a particular node i and to determine the 
probability of have K randomly chosen nodes linked to 
it. The probability that a link comes from set £{t) with 
the corresponding node having value 1 is z\ = E\/N . If 
the link comes from the same set but the node of £ (t) has 
value 0, the probability will be zq = Eo/N. wq = Fq/N 
and wi = F\/N are the probabilities that a link comes 
from F(t) when the corresponding elements of the set 
have values and 1, respectively. Since E + E\ + F + 
Fi = N, it is obvious that zq + z\ + wq + w\ = 1. 

We are interested in the e valuation of W\ , the proba- 
bility that Oi(t + 1) = 1 and a -h (t + 1) = 0, and of Wq, 
the probability that criit + 1) = and Oi(t + 1) = 1. The 
first step is to study the situation in which the node i has 



K — 1 links in set £{t) and only one link in ^(t). The list 
of the probabilities of the possible link configurations is 



p(l) 



K 



1 

K -1 
2 



wo, 



K 



,K-1 I K 1 \ K-2, 



Z K-3 Z 2 
z z li 



~K-l 



(10) 



where we must multiply each element of the first list by 
each element of the second one. We must now evalu- 
ate the output of each possible configuration under the 
operator XOR. For the first configuration 
the XOR operation furnishes that <jj = 1 in configura- 
tion {a(t + 1)} and ai — in {a(t + 1)}, implying that 
p(^) wiz^ -1 will contribute to W\. Note that the extra 
probability p is related to the XOR operator. The sec- 
ond configuration ((^) wi z$ ~ 2 Z\) will contribute 
to Wq because the application of XOR g ive- us that cr; = 
in {a(t + 1)} and Ui = 1 in {a(t + 1)}. Since the third 
term ((f)u>i ^^z^^zf) will contribute to Wi, it is 
easy to infer that for configurations beginning with w\, 
the terms z™ with m even contribute to W\ and the ones 
with to odd enter in Wq. A similar analysis shows that 
for configurations beginning with wq, the terms with m 
even contribute to Wq and the ones with to odd enter in 
W\. When we apply the AND operator, only the terms 
wiz^ 1 and wqz^~ x give no null contributions to W\ 
and Wq, respectively. Therefore the contributions of the 
list (fTU)) to W\ and Wq can be written as 



pw\ 



K-l 

E 

m=0 
m even 



K - 1 



K — rn. m 



m=l 
m odd 



rn 



w^ = 



PWQ 



K-l 

E 

m even 



K - 1 



^ (K V 



m— 1 
m odd 



. K — rn m 



l zT + qwQzi 



Using Equations ([6]) and (O, w[ 1%> can be written as 



(i) 



-(Wl + Wq)(zq + Zi) 



K-l 



(wi - wq)(z - zi) +qwiz 1 



The equation for Wq is obtained by changing wq by w\, 
and Wi by Wq in the above equation. 

The second step is to study the situation in which the 
node i has K — 2 links in set £ (t) and two links in F{t). 



.K-i 
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Now, the list of the probabilities of the possible link con- 
figurations is 



Pf 3 = (f ){wl 2w 0Wl , W \}{4 



K-1 



(11) 



K - 2 



K - 2 



z K-4 2 



■ ; 2 1 



There is no contribution of the XOR operator. The AND 
operator furnishes that only two configurations {w\z^ ~ 2 



and 



n 2 z K ~ 2 
u z 1 



), multiplied by probability q, contribute to 



Wi and Wq. Therefore we find that 



(2) 



K 2 \wlz«- 2 , and WP 



K\ _ 
2 ) qW " ' [ 



2 K-1 



The third step is to evaluate the probabilities gener- 
ated by the application of XOR and AND in the situation 
in which the node i has K — 3 links in set £{t) and three 
links in J-{t). This situation is similar to the first one. 
We obtain that 



p, 



|(wi+W ) 3 (^0 + 2l) K 3 



3_A-3 



+ |K -WoY{zq-Zi)* i +qwlz{ 



The result for W is identical with the previous one if 
we change wq by Wi, and w\ by wq. 

The fourth step is similar to the second one, and so on. 



Since W\ — J2m=i Wi™' , we have that 



(m) 



^ = 1 £ 



A' 



odd 



K 



[(wi +w ) m (z + z 1 ) 



K—m 



A 



+ {w 1 - w Q ) m (z - zi ) 



K — rn 1 



m — 1 



1 



To obtain Wo we substitute wq by wi, and wi by wq in 
this equation. Using again Eq. ©, we obtain that 



Wi 



q[( Wl + Zl ) K - zf ] + |{1 - [1 - 2(t«i + U*)] 



A' 



+ [l-2(z 1 + Wo )] K -[l-2(z 1 + Wl )] K }, (12) 



W = q[(w + z x ) K - z* } + |{1 - [1 - 2( Wl + w )] K 
+ [l-2{z 1 +w 1 )] K - [l-2( Zl +w )] K }. (13) 



Observe that the equations for W\, Wo and Z\ describe 
completely our system, since the equation for Zo is ob- 
tained from the normalization condition. However it is 
usual to work with variables M(p,t + 1), the density of 
l's and with D(p,t + 1), the Hamming distance. From 
the definition of the Hamming distance we have that 
D(p,t + 1) = wi, t +i + w ot +i, implying that 



D(p,t + 1) = q{M(p,t) K +[D(p,t)-M(p,t) 
- 2z^ t } + ^{l-[l-2D(p,t)} K }. 



2zi 



lA 



The equations for D(p, t + 1), M(p, t + l) and zi : t+i also 
describe completely the system. However, if w± = wq 
we can see from Eqs. JT2J! and (JTSJ) that W\ = W . 
Solving numerically these equations, we obtain that each 
initial configuration with w\ ^ wq evolves to a fixed point 
with W\ = Wo. Then we can assume that wq = wi, 
without loss of generality, and the dynamics of the system 
is described by only two equations, namely 



2q{M(p,t)K-[M(p,t) 



D(p,t+1) 
M(p,t+l) = qM(p,t) 



{l-(l-2D(p,t)) K } 



A 



(14) 
(15) 



£{l-[l-2D(p,t)] K }. 



Finally let us consider the model with the scale-free 
topology. Since each node now have k links with prob- 
ability P(k), the above equations, which are valid for 
w\ = wo, can be written as 



D(jp,t + 1) 




M(p,t + 1) 



'~2~ 1 J 

2 ^P(fc)[l-2 J D(p,t)] fc , (16) 

k=1 



l-^P(fc)[l-2Af(p,i)] A 



oc 

g^F(fc)Af(p,O fc . 

k=1 



(17) 



Assuming homogeneity, we can identify Wi as W\ t t+\i the 
fraction of l's of set Fit) at time t + l, and Wo as wo.t+i- 
Note that the fraction of l's of the system is given by 
M(p,t + 1) = zi.t+i + wi,t+i and that it was already 
evaluated (see Eq.®). Identifying z\t+i with Z\, we 
obtain that 

Z x = M(p,t + l) - Wi. 



C. M{p,t) and D(jp,t) for k mm = 1 

The main difference between this case and the previ- 
ous one is the self-regulation mechanism: the node itself 
always participates in its own dynamics. Itself and the K 
nodes connected to it are the control elements of the dy- 
namics. Following a similar procedure of the subsection 
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IIV A[ we obtain that 



M(p,t + 1) = | I 1 - j^P(k)[l - 2M(p,t)] k+1 



fe=i 



(18) 



fc=i 



where P(k) is the probability that a node had fc links. 

The evaluation of the Hamming distance follows sim- 
ilar steps of subsection IIVBI For scale-free systems we 
again have the quantities Wo, W\ and Z\. It turns out 
that when w% = Wq we have W\ = Wq. When w\ — wq, 
the dynamics of the system is given by Eq. (TTS)) and by 

( oo 

D(p,t+l) = 2q\Y J P(k){M{p,t) k+1 



,fc=i 



[M(p,t) 



D(jp,t) 



fe+i 



+ f -f £i'(A)[l-2i?(p,t)]* +1 .(19) 
fc=i 

If we put P(fc) = in Eqs. (TT8"|) and ([K))) . we obtain 
the results for the model without a scale-free topology. 
They are similar to the ones obtained for the case with- 
out self-regulation, but with K replaced by K + 1 (see 
Eqs. (|14|) and (fT5|)). For scale-free systems, the dynam- 
ics with self-regulation is similar to the usual dynamics 
if we change k by k + 1, except in the distribution of 
connectivity P(k) (see Eqs. ([TBI and ITT]) . Therefore the 
dynamics with self-regulation is different from the usual 
case. Let us investigate if this fact is also true for the 
Kauffman NK model. 



probability of having different values in {a(t + 1)} and 
{a(t + 1)}. Due to the random Boolean function assign- 
ment any node can be or 1 with probability 1/2, and 
the probability that a^t + l) will have different values in 
{a(t+l)} and {a(t + 1)} is 1/2. Therefore the probabil- 
ity W is given by 



W 



z K - x w 



z K - 2 w 2 



w 



K 



Assuming that the system is homogeneous, we can iden- 
tify the fraction of nodes with different values in {a(t+l)} 
and {a(t + 1)}, W, with the Hamming distance D(p,t + 
1). Using that w = D(p,t) in the previous equations we 
obtain the traditional equation of Derrida and Pomeau 
[30| . namely 

D(p,t + l) = ^{l-[l-D(p,t)] K }. 

Let us consider now a model with self-regulation. 
Moreover, each node also has K links connected to any 
of the N nodes of the network. We focus on node i. This 
node has probabilities = z and w% = w to belong to 
sets £{t) and .F(i), respectively. We want again to com- 
pute W. If node i is in set T(t) it has probability 1/2 
to contribute to W, independently of the K links. Oth- 
erwise, at least one of the K links must be in set F{t). 
Taken in account these two situations, we have that 



Wj Zt 

W = Y + 1 



+ 10 



A" 



Identifying W with D(p,t + 1) and w with D(p,t), we 
find that the Hamming distance is given by 



1 



D(p,t + 1) = -{1- [1 _D(p, <)]*+*} 



(20) 



D. Kauffman model with self-regulation 

The Kauffman NK model consists of N nodes holding 
logical variables Ui. Each node is connected with any K 
nodes of the network. Observe that a node i can have a 
link to itself with small probability (K/N) . The dynam- 
ics, given by Eq. (JlJ, is determined by a random Boolean 
function Gi. In the Derrida and Pomeau annealed ap- 
proximation [3(3], the Boolean function and the K nodes 
are randomly chosen at each time step. The configuration 
{&(£)} is split in the sets T(t), which consists of nodes 
having different values of a in configurations {cr(t)} and 
{a(t)}, and £{t) when the previous condition does not 
hold. Then we are able to define the probabilities w and 
z that a link of a particular node comes from sets 
and £ (t), respectively. Obviously we have that z + w = 1. 
We want to evaluate the probability W that nod e i will 
have different values in {a(t + 1)} and {a(t + 1)}. If all 
K links came from £ (i), erj(i + 1) will have the same 
value in both {a(t + 1)} and {a(t + 1)}. However, if at 
least one link comes from J-(t), (Ti(t + 1) has a positive 



Observe that again this expression is similar to the pre- 
vious one if K is replaced by K + 1 . 

Both models can be studied in a scale-free topology. It 
is easy to obtain that 



D(p,t + \) 



i|l-^P(fc)[l-^,t)] e |, (21) 



where 8 = k + 1 for the case with self-regulation, and 
6 = k otherwise. These results are easily generalized to 
taken in account the Derrida parameter pd (see Derrida 
and Pomeau [3(|). In this case the probability 1/2 of 
Eqs. (|20|) and (j2Tj) must be replaced by the corresponding- 
probability 2prf(l - pa). 



V. MEAN-FIELD AND SIMULATION RESULTS 

A. Results for k m in ~ 2 

The fixed points for the model without a scale-free 
topology are obtained by putting M(p, t+1) = M (p, t) = 
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TABLE I: Values of density of l's (M(p)) and Hamming dis- 
tance (D(p)) for the model without a scale- free topology with 
kmin — 2, N — 10000, and the dynamics described by XOR 
and AND functions without self-regulation. The subscript 
sim is the simulated result and the ann refers to annealed 
(MF) solution. The error in the last digit of an evaluated 
quantity in the simulations is between parentheses. 



K 


P 








■D sim 


2 


0.0 


0.000 


000(11 


0.000 


000(21 




0.2 


0.000 


OOOfD 


0.000 


000(2) 




0.5 


0.002 


0.054(3) 


0.002 


0.060(3) 




0.7 


0.364 


0.362(3) 


0.399 


0.397(2) 




1.0 


0.500 


0.500(2) 


0.500 


0.500(1) 


5 


0.0 


0.000 


0.000(1) 


0.000 


0.000(1) 




0.2 


0.000 


0.017(3) 


0.000 


0.016(2) 




0.5 


0.241 


0.241(2) 


0.242 


0.241(2) 




0.8 


0.402 


0.402(2) 


0.404 


0.404(2) 




1.0 


0.500 


0.500(1) 


0.500 


0.500(2) 


10 


0.0 


0.000 


0.000(1) 


0.000 


0.000(1) 




0.2 


0.084 


0.084(1) 


0.084 


0.083(2) 




0.5 


0.250 


0.249(1) 


0.250 


0.250(2) 




0.8 


0.400 


0.400(1) 


0.400 


0.400(1) 




1.0 


0.500 


0.500(1) 


0.500 


0.500(1) 



M* in the map given by Eq. ([5]) . The fixed point M* = 
always exists, and the local stability parameter A, given 
by 

A = dMjpt + 1) = K _ x + km k-x 

dM[p,t) Mr 

tell us that M* = is stable for A = pK < 1. It means 
that lim{->oo M(p, t) — for any initial value M(p, 0). 
When pK > 1 the initial conditions are attracted to a 
non null fixed point. It means that in the (p, K) plane, 
there is a curve, given by equation pK = 1, separating 
the region in which M* = is stable from the one that 
M* = is not stable. These features can be easily illus- 
trated for K = 2. In this case, the non null fixed point is 
given by M* = (1 — 2p)/(l — 3p) and the local stability 
parameter, evaluated for non null M», is A = 2(1 — p). 
Then we have that A < 1 for p > 1/2. It implies that 
for p > 1/2, the non null fixed point is attractive. For 
K > 2 the evaluation of the fixed point and of A were per- 
formed numerically. In table (JH we compare the results 
of the annealed approximation for three values K with 
the ones obtained by numerical simulations of N = 10 4 
nodes, with the average quantities evaluated after 10 3 
mcs in 3000 samples. In the simulation results, the num- 
bers between parentheses are the errors that affect the 
last digits. Similar results were obtained for other values 
of p. We can conclude that the MF results for the frac- 
tion of l's agree very well with those from the numerical 
simulations. 



TABLE II: Values of density of l's (M(p)) and Hamming dis- 
tance (D(p)) for Boolean dynamics described by XOR and 
AND functions without self-regulation on scale-free networks 
with N = 10000, k min = 2, and r = 0.5 (7 = 3). The sub- 
scripts sim and ann refer to the simulated and MF results, 
respectively. The errors evaluated in the simulations are be- 
tween parentheses. 



p 










0.2 


0.003 


0.038(3) 


0.038 


0.007(3) 


0.3 


0.075 


0.085(3) 


0.113 


0.039(3) 


0.5 


0.227 


0.216(3) 


0.242 


0.178(3) 


0.7 


0.352 


0.350(4) 


0.353 


0.340(4) 


0.8 


0.406 


0.409(4) 


0.404 


0.414(4) 



The fraction of l's in a scale-free network is described 
by Eq. ([§]). Again the fixed point M* = is always 
present. A can also evaluated and we have that the M* = 
is stable if p < l/(k). When M* = is not stable, we 
obtain numerically that there is a non null fixed point 
attracting all the initial conditions. These regions are 
separated in the (p, (fc)) plane by a curve described by 

p(k) = 1. (22) 

In Tab. (JTTJ) we compare the results of the annealed ap- 
proximation with the ones obtained from numerical sim- 
ulations with 10 4 nodes and r — 0.5 (7 = 3). The M* 
results obtained by MF solution agree well with the ones 
from simulations. Similar results are obtained for other 
values of r. 

The Eqs. (fl"4| and (|16|) for the Hamming distance can 
be numerically solved to furnish the fixed points in the 
cases of the models without and with scale-free topology. 
Note that we are using that wi = Wq in both cases. In 
Tabs. (P) and (JTl]) we can compare the results obtained 
from the annealed approximation with those obtained by 
the numerical simulations. We see that they agree well. 
This conclusion holds for other values of the parameter 

P- 

B. Results for fc m i„ = 1 

In this case we have self-regulation. The fixed points 
for M* and D„, obtained from Eqs. (| and (|T5)) with 
P(k) — 5k, k, are displayed in Tab. (|III|) . We can see 
that MF results are different from the ones obtained from 
numerical simulations for small K (K = 1 and K = 3), 
although they are similar when K is large (K = 10). 
In order to check our analytical approximation, we have 
also performed numerical simulations with an annealed 
dynamics. At each mcs we have randomly chosen the K 
nodes connected with each node of the network. As we 
can see, the numerical results agree very well with the 
ones obtained from MF. 
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TABLE III: Values of density of l's (M(p)) and Hamming 
distance (D(p)) for the model without a scale- free topology 
with k m i n = 1, AT = 10000, and the dynamics described by 
XOR and AND functions with self-regulation. The subscript 
aim is the simulated result and the ann refers to annealed 
(MF) solution. 



! \ 


P 


M 


M ■ 


n 


n 

i-* srrn 


1 


0.0 


0.000 


0.002(2) 


0.000 


1 1 /will / 1 \ 

0.000(3) 




0.2 


0.000 


0.130(3) 


0.000 


0.058(3) 




0.5 


0.002 


0.277(4) 


0.002 


0.114(4) 




0.8 


0.429 


0.406(5) 


0.454 


0.172(5) 




1.0 


0.500 


0.497(4) 


0.500 


0.382(5) 


3 


0.0 


0.000 


0.000(2) 


0.000 


0.000(3) 




0.2 


0.000 


0.100(4) 


0.000 


0.073(4) 




0.5 


0.230 


0.250(3) 


0.232 


0.233(3) 




0.8 


0.405 


0.400(3) 


0.410 


0.399(3) 




1.0 


0.500 


0.500(3) 


0.500 


0.500(3) 


10 


0.0 


0.000 


0.000(2) 


0.000 


0.000(3) 




0.2 


0.088 


0.010(2) 


0.088 


0.099(3) 




0.5 


0.250 


0.250(3) 


0.250 


0.250(2) 




0.8 


0.400 


0.400(3) 


0.400 


0.400(3) 




1.0 


0.500 


0.500(3) 


0.500 


0.500(3) 



The same conclusions hold for the scale-free topology. 
By using Eq. (|18l) in the analysis of stability of the Af* = 
fixed point, we obtain that the curve separating the two 
regions is given by 

p«fc) + l) = l. (23) 

Even if we take into account that (k) of the above equa- 
tion begins with k = 1 and the one of Eq. (|22|) begins 
with k = 2, Eqs. (22|) and ([23]) are different. This implies 
that self-regulation changes the dynamical behavior. 

Table ITVl shows the results of computational simulation 
and the annealed approximation. It can be seen that the 
values of M and D obtained via MF when fc m , n = 1 are a 
bit smaller than the ones of the k m i n = 2 case. Another 
important feature is the relative good agreement between 
the values of M obtained via simulation and that from 
annealed approximation. This match does not occur for 
the Hamming distance. As we can see in Tab. IIV1 the 
spreading damage calculated via MF is quite bigger than 
the one obtained in simulations. We can conclude by 
analyzing the Hamming distance, that the self-regulated 
nodes introduce a new dynamical behavior, with distinct 
properties when compared to the behavior of non-self- 
rcgulated ones. Is worth to comment that self- regulation 
is a common feature in biological networks. Maybe it is 



a process used in order to increase homeostasis, reducing 
the effect of a damage introduced in the system. 

TABLE IV: Values of density of l's (M(p)) and Hamming 
distance (D(p)) for Boolean dynamics described by XOR and 
AND functions with self-regulation on scale-free networks 
with N = 10000, k min = 1, and r = 0.5 (7 = 3). The 
subscripts sim and ann refer to the simulated and MF re- 
sults, respectively. The errors evaluated in the simulations 
are between parentheses. 



V 








-Dsim 


0.2 


0.091 


0.097(3) 


0.090 


0.000(3) 


0.3 


0.141 


0.140(2) 


0.140 


0.000(3) 


0.5 


0.246 


0.228(3) 


0.243 


0.001(3) 


0.7 


0.351 


0.325(2) 


0.345 


0.013(4) 


0.8 


0.401 


0.380(2) 


0.396 


0.042(4) 



VI. SUMMARY 

In this work we studied Boolean dynamics in Kauff- 
man models and in scale-free networks. The dynamical 
models assigned only XOR and AND operators to the 
nodes with probability p and 1 — p. Regarding the in- 
puts of the above cited Boolean networks, two types of 
dynamics were used. In the first one, the state of the 
nodes was regulated by the state of all nodes connected 
to them. The second type was similar to the first one, 
with the difference that the state of the node was used as 
its own input. Thus, in the first case we did not have self- 
regulation as in the second one. As shown in the results, 
these two types of dynamics presented quite different be- 
haviors. In both cases a computational simulation and 
an analytical mean-field approximation were performed 
in order to compare the density of l's, namely M, and 
the Hamming distance D. The results for the dynam- 
ics with no self-regulation generated good agreement be- 
tween simulations and the MF approach. However, the 
case with self-regulation had a clear disagreement with 
respect to D and M, for small values of p. 



The authors thank J. F. F. Mendes for useful discussions. 
We thank the referees for useful suggestions. ACS and 
JKLS thank to Fundacao de Amparo a Pesquisa de MG 
(FAPEMIG), a Brazilian agency, for partial financial 
support. JKLS thanks to Conselho Nacional de Pesquisa 
(CNPq) for partial financial support. 
*Electronic address: alcidescs@gmail.com 
^Electronic address: jaff@fisica.ufmg.br 



[1] A. L. Barabasi and R. Albert, Science 286, 509 (1999). [2] A. L. Barabasi and Z. N. Oltvai, Nature Rev. Genet. 5, 



101 (2004). 

[3] S. N. Dorogovtsev and J. F. F. Mendes, Adv. Phys. 51, 
1079 (2002). 

[4] L. Wardil and J. K. L. da Silva, Braz. J. Phys. 38, 350 
(2008). 

[5] A. T. Bernardes, D. Stauffer and J. Kertesz, Eur. Phys. 

J. B 25, 123 (2002). 
[6] D. J. Watts and S. H. Strogatz, Nature 393, 440 (1998). 
[7] S. Redner, Eur. Phys. J. B 4, 131 (1998). 
[8] B. A. Huberman, P. L. T. Pirolli, J. E. Pitkow and R. 

M. Lukose, Science 280, 95 (1998). 
[9] R. Albert, H. Jeong and A. L. Barabasi, Nature 401, 130 

(1999). 

[10] B. A. Huberman and L. A. Adamic, Nature 401, 131 
(1999). 

[11] M. Barthelemy and L. A. N. Amaral, Phys. Rev. Lett. 

82, 5180 (1999). 
[12] M. E. J. Newman and D. J. Watts, Phys. Lett. A 263, 

341 (1999). 

[13] A. Barrat and M. Weigt, Eur. Phys. J. B 13, 547 (2000). 
[14] L. A. N. Amaral, A. Scala, M. Barthelemy and H. E. 

Stanley, Proc. Natl. Acad. Sci. USA 97, 11149 (2000). 
[15] L. A. Barbosa, A. Castro e Silva and J. K. L. da Silva, 

Phys. Rev. E 73, 41903 (2006). 
[16] H. Jeong, B. Tombor, R. Albert, Z. N. Oltvai and A. L. 



11 



Barabasi, Nature 407, 651 (2000). 
[17] S. A. Kauffman, J. Theor. Biol. 22, 437 (1969). 
[18] A. Castro e Silva, J. K. L. da Silva and J. F. F. Mendes, 

Phys. Rev. E 70, 66140 (2004). 
[19] K. Iguchi, S. I. Kinoshita and H. S. Yamada, J. Theor. 

Biol. 247, 138 (2007). 
[20] E. Novikov and E. Barillot, BMC Sis. Biol. 2, 8 (2008). 
[21] B. Drossel and F. Greil, Phys. Rev. E 80, 026102 (2009). 
[22] M. Aldana, Physica D 185, 45 (2003). 
[23] D. S. Lee and H. Rieger, J. Phys. A: Math. Theor. 41, 

415001 (2008). 

[24] H. J. Zhou and R. Lipowsky, Proc. Natl. Acad. Sci. USA 
102, 10052 (2005). 

[25] S. Kauffman At Home in the Universe, Oxford University 
Press, New York, 1995. 

[26] I. Shmulevich, E. R. Dougherty and W. Zhang, Proc. 
IEEE 90, 1778 (2002). 

[27] N. Guelzim, S. Bottani, P. Bourgine and F. Kepes, Na- 
ture Genet. 31, 60 (2002). 

[28] P. L. Krapivsky and S. Redner, Phys. Rev. E 63, 066123 
(2001). 

[29] F. Greil and B. Drossel, Eur. Phys. J. B 57, 109 (2007). 
[30] B. Derrida and Y. Pomeau, Europhys. Lett. 1, 45 (1986). 



